/*
 * Copyright (c) 2009-2013, Peter Abeles. All Rights Reserved.
 *
 * This file is part of Efficient Java Matrix Library (EJML).
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package mikera.matrixx.decompose.impl.svd;

import mikera.matrixx.AMatrix;
import mikera.matrixx.Matrix;
import mikera.matrixx.decompose.Bidiagonal;
import mikera.matrixx.decompose.IBidiagonalResult;
import mikera.vectorz.AVector;
import mikera.vectorz.Vector;

/**
 * <p>
 * Computes the Singular value decomposition of a matrix using the implicit QR algorithm
 * for singular value decomposition.  It works by first by transforming the matrix
 * to a bidiagonal A=U*B*V<sup>T</sup> form, then it implicitly computing the eigenvalues of the B<sup>T</sup>B matrix,
 * which are the same as the singular values in the original A matrix.
 * </p>
 *
 * <p>
 * Based off of the description provided in:<br>
 * <br>
 * David S. Watkins, "Fundamentals of Matrix Computations," Second Edition. Page 404-411
 * </p>
 *
 * @author Peter Abeles
 */
public class SvdImplicitQr {

    private int numRows;
    private int numCols;

    // dimensions of transposed matrix
    private int numRowsT;
    private int numColsT;

    // if true then it can use the special Bidiagonal decomposition
//    private boolean canUseTallBidiagonal;

    // If U is not being computed and the input matrix is 'tall' then a special bidiagonal decomposition
    // can be used which is faster.
    private IBidiagonalResult bidiagResult;
    private SvdImplicitQrAlgorithm qralg = new SvdImplicitQrAlgorithm();

    double diag[];
    double off[];

    private Matrix Ut;
    private Matrix Vt;

    private double singularValues[];
    private int numSingular;

    // compute a compact SVD
    private boolean compact;
    // What is actually computed
//    private boolean computeU;
//    private boolean computeV;

    // What the user requested to be computed
    // If the transpose is computed instead then what is actually computed is swapped
//    private boolean prefComputeU;
//    private boolean prefComputeV;

    // Should it compute the transpose instead
    private boolean transposed;

    // Either a copy of the input matrix or a copy of it transposed
    private Matrix A_mod = Matrix.create(1,1);
    
    public static SVDResult decompose(AMatrix A, boolean compact) {
    	SvdImplicitQr svd = new SvdImplicitQr(compact);
    	return svd._decompose(A);
    }

    /**
     * Configures the class
     *
     * @param compact Compute a compact SVD
     * @param computeU If true it will compute the U matrix
     * @param computeV If true it will compute the V matrix
     */
    /*Once BidiagonalDecomposition is implemented, use the commented out constructor
    and add:
    "@param canUseTallBidiagonal If true then it can choose to use a tall Bidiagonal decomposition to improve runtime performance."
    to doc*/
//    public SvdImplicitQr(boolean compact, boolean computeU, boolean computeV,
//    		boolean canUseTallBidiagonal )    
    public SvdImplicitQr(boolean compact) {
        this.compact = compact;
//        this.prefComputeU = computeU;
//        this.prefComputeV = computeV;
//        this.canUseTallBidiagonal = canUseTallBidiagonal;
    }

    public AVector getSingularValues() {
        return Vector.wrap(singularValues);
    }

    public int numberOfSingularValues() {
        return numSingular;
    }

    public boolean isCompact() {
        return compact;
    }

    public AMatrix getU() {
//        if( !prefComputeU )
//            throw new IllegalArgumentException("As requested U was not computed.");
    	return Ut.getTranspose();
    }

    public AMatrix getV() {
//        if( !prefComputeV )
//            throw new IllegalArgumentException("As requested V was not computed.");
        	return Vt.getTranspose();
    }

    public AMatrix getS() {
        int m = compact ? numSingular : numRows;
        int n = compact ? numSingular : numCols;

        Matrix S = Matrix.create(m,n);
        
        for( int i = 0; i < numSingular; i++ ) {
            S.unsafeSet(i,i, singularValues[i]);
        }

        return S;
    }

    public SVDResult _decompose(AMatrix _orig) {
//    	Creating a copy so that original matrix is not modified
    	Matrix orig = _orig.copy().toMatrix();
        setup(orig);

        if (bidiagonalization(orig)) {
            return null;
        }

        if( computeUSV() )
            return null;

        // make sure all the singular values or positive
        makeSingularPositive();

        // if transposed undo the transposition
        undoTranspose();

        return new SVDResult(getU(), getS(), getV(), getSingularValues());
    }

    private boolean bidiagonalization(Matrix orig) {
        // change the matrix to bidiagonal form
        if( transposed ) {
            A_mod = orig.getTransposeCopy().toMatrix();
        } else {
            A_mod = orig.copy().toMatrix();
        }
        bidiagResult = Bidiagonal.decompose(A_mod, compact);
        return bidiagResult == null;
    }

    /**
     * If the transpose was computed instead do some additional computations
     */
    private void undoTranspose() {
        if( transposed ) {
            Matrix temp = Vt;
            Vt = Ut;
            Ut = temp;
        }
    }

    /**
     * Compute singular values and U and V at the same time
     */
    private boolean computeUSV() {
        diag = bidiagResult.getB().getBand(0).toDoubleArray();
        off = bidiagResult.getB().getBand(1).toDoubleArray();
        qralg.setMatrix(numRowsT,numColsT,diag,off);

//        long pointA = System.currentTimeMillis();
        // compute U and V matrices
//        if( computeU )
            Ut = bidiagResult.getU().getTranspose().toMatrix();
//        if( computeV )
            Vt = bidiagResult.getV().getTranspose().toMatrix();

        qralg.setFastValues(false);
//        if( computeU )
            qralg.setUt(Ut);
//        else
//            qralg.setUt(null);
//        if( computeV )
            qralg.setVt(Vt);
//        else
//            qralg.setVt(null);

//        long pointB = System.currentTimeMillis();

        boolean ret = !qralg.process();

//        long pointC = System.currentTimeMillis();
//        System.out.println("  compute UV "+(pointB-pointA)+"  QR = "+(pointC-pointB));

        return ret;
    }

    private void setup(Matrix orig) {
        transposed = orig.columnCount() > orig.rowCount();

        // flag what should be computed and what should not be computed
        if( transposed ) {
//            computeU = prefComputeV;
//            computeV = prefComputeU;
            numRowsT = orig.columnCount();
            numColsT = orig.rowCount();
        } else {
//            computeU = prefComputeU;
//            computeV = prefComputeV;
            numRowsT = orig.rowCount();
            numColsT = orig.columnCount();
        }

        numRows = orig.rowCount();
        numCols = orig.columnCount();

        diag = new double[ numColsT ];
        off = new double[ numColsT-1 ];

        // if it is a tall matrix and U is not needed then there is faster decomposition algorithm
//        if( canUseTallBidiagonal && numRows > numCols * 2 && !computeU ) {
//            if( bidiag == null || !(bidiag instanceof BidiagonalDecompositionTall) ) {
//                bidiag = new BidiagonalDecompositionTall();
//            }
//        } else if( bidiag == null || !(bidiag instanceof BidiagonalDecompositionRow) ) {
//            bidiag = new BidiagonalDecompositionRow();
//        }
//        TODO: ^Choose between BidiagonalTall and BidiagonalRow once  BidiagonalTall
//        is implemented
    }

    /**
     * With the QR algorithm it is possible for the found singular values to be negative.  This
     * makes them all positive by multiplying it by a diagonal matrix that has
     */
    private void makeSingularPositive() {
        numSingular = qralg.getNumberOfSingularValues();
        singularValues = qralg.getSingularValues();

        double[] UtData = Ut.asDoubleArray();
        for( int i = 0; i < numSingular; i++ ) {
            double val = qralg.getSingularValue(i);

            if( val < 0 ) {
                singularValues[i] = 0.0d - val;

//                if( computeU ) {
                    // compute the results of multiplying it by an element of -1 at this location in
                    // a diagonal matrix.
                    int start = i* Ut.columnCount();
                    int stop = start+ Ut.columnCount();

                    for( int j = start; j < stop; j++ ) {
                        UtData[j] = 0.0d - UtData[j];
//                    }
                }
            } else {
                singularValues[i] = val;
            }
        }
    }
    
    public int numRows() {
        return numRows;
    }

    public int numCols() {
        return numCols;
    }
}
